Modeltime with tidymodels

tidymodels를 이용한 시계열 모델링

Introduction

modeltime 패키지는 tidymodels와 연동이 가능한 시계열 모델링 관련 패키지이다. tidymodels에도 auto arima, ma 모형 같은 간단한 시계열 모델이 있지만 다른 머신러닝 모델에 비해서 쓸 수 있는 모델이 제한적이다. 이러한 단점을 modeltime 패키지가 해결해준다.

modeltime 시계열 모델링에 특화된 패키지로 시계열에 특화된 모델링과 전처리 관련 함수가 내장되어 있고, tidymodels의 워크플로우를 거의 그대로 이용할 수 있어서 향후 tidymodels가 R의 대표적인 머신러닝 패키지가 된다면 시계열 파트에서는 modeltime 패키지가 주축이 될 것 같다.

modeltime 패키지는 하나의 단일 패키지가 아니라 다양한 머신러닝 패키지와 연동해서 하나의 시계열 생태계를 구축하고 있다. 대표적인 패키지는 다음과 같다.

• Modeltime : 시계열 머신러닝 관련 패키지

• Modeltime H2O : H2O의 autoML과 연동이 가능함

• Modeltime GluonTS : 시계열 관련 딥러닝 패키지

• Modeltime Ensemble : Modeltime 관련 앙상블 패키지

• Modeltime Resample : Backtesting 관련 패키지

• Timetk : feature engineering, Data wrangling, time series visualization

Modeltiem workflow

Modeltime 패키지는 시계열 모델링에 맞는 고유한 워크플로우가 존재하며, tidymodels와 마찬가지로 이 워크플로우 그대로 진행을 해야만 에러 없이 동작한다. 사실 실제 모델링할 때와 동일한 매커니즘이어서 어색함 없이 따라할 수 있다.

workflow의 순서는 다음과 같다.

1. create modeltime table

   • 다양한 모델을 training data에 fitting하고 modeltime table에 저장하는 단계

   • 모델이 재대로 fitting 되었는지 확인 및 이후 워크플로우에 맞는 table 구성

2. calibrate

   • modeltime table에 저장된 모델을 test 데이터에 fitting하는 단계

   • test data에 fitting함으로써 모델의 성능 파악 가능

3. Refit

   • training data, test data를 합친 original data에 refit하는 단계

   • refit한 모델을 이용해서 forecast를 진행

Model

Modeltime 패키지는 다양한 시계열 모델을 지원하는데 대표적인 모델은 다음과 같다.

• auto_arima_boost: a automatic ARIMA boosted model specification
• arima_reg : arima regression model
• exp_smoothing : automated exponential smoothing
• prophet_boost: a prophet model specification
• prophet_boost_log: a prophet model specification with logarithmic growth
• mars: a Multivariate Adaptive Regression Splines model specification
• nnetar: a Neural network time series model specification

Modeltime 패키지의 한 가지 단점은 패키지를 만든 저자가 business-science라는 course를 운영 중인 것 같은데 modeltime 관련 세부 자료는 유료 course에서 공개를 하는 것 같다… 하지만 R CRAN에도 등록이 되어있고 튜토리얼도 제공되기 때문에 앞으로 사용하는 사람이 많아지면 관련 자료도 많아질 것 같다.

Modeltime tutorial

데이터는 데이콘에서 제공하는 발전량 데이터 일부를 뽑아서 활용했다.

Preparations (준비작업)

Libraries

library(tidyverse)
library(tidymodels)
library(lubridate)
library(data.table)
library(skimr)
library(timetk)
library(modeltime)
library(gt)
library(timetk)
library(tidyquant)
library(sknifedatar)

# Visualization
library(ggthemes)
library(ggsci)
library(viridis)
library(ggExtra)


theme_set(theme_bw())

Data load

file_path <- getwd()
files <- list.files(file_path)
files

[1] "modeltime-with-tidymodels.html" 
[2] "modeltime-with-tidymodels.Rmd"  
[3] "modeltime-with-tidymodels_files"
[4] "rdata.csv"                      

rdata <- fread(file.path(file_path, "rdata.csv"))

Data overview (데이터 기본정보)

Data

head(rdata)

                  time 기온 강수 풍속 풍향 습도 증기압 이슬점   기압
1: 2015-01-01 01:00:00 -4.4 0.00  5.4  340   47    2.1  -14.0 1020.3
2: 2015-01-01 02:00:00 -4.6 0.00  4.9  340   50    2.2  -13.4 1020.3
3: 2015-01-01 03:00:00 -4.7 0.05  6.2  320   50    2.2  -13.5 1020.7
4: 2015-01-01 04:00:00 -5.0 0.00  5.0  320   56    2.4  -12.4 1020.6
5: 2015-01-01 05:00:00 -5.0 0.00  5.5  320   52    2.2  -13.3 1020.4
6: 2015-01-01 06:00:00 -5.3 0.00  4.0  340   58    2.4  -12.2 1020.7
   해면기압 일조 일사 전운량 시정 지면온도 floating warehouse dangjin
1:   1024.0   NA    0      6 1500     -4.4       NA         0       0
2:   1024.0   NA    0     NA 1750     -4.6       NA         0       0
3:   1024.5   NA    0      6 2000     -4.7       NA         0       0
4:   1024.4   NA    0      6 2000     -5.0       NA         0       0
5:   1024.2   NA    0      6 2000     -5.0       NA         0       0
6:   1024.5   NA    0      6 2000     -5.3       NA         0       0
   ulsan hour
1:     0    1
2:     0    2
3:     0    3
4:     0    4
5:     0    5
6:     0    6

glimpse(rdata)

Rows: 55,392
Columns: 20
$ time      <dttm> 2015-01-01 01:00:00, 2015-01-01 02:00:00, 2015-01~
$ 기온      <dbl> -4.4, -4.6, -4.7, -5.0, -5.0, -5.3, -5.7, -5.7, -5.4~
$ 강수      <dbl> 0.00, 0.00, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00~
$ 풍속      <dbl> 5.4, 4.9, 6.2, 5.0, 5.5, 4.0, 5.0, 4.6, 6.7, 6.5, 6.~
$ 풍향      <dbl> 340, 340, 320, 320, 320, 340, 340, 340, 320, 320, 32~
$ 습도      <dbl> 47, 50, 50, 56, 52, 58, 58, 56, 57, 59, 60, 62, 58, ~
$ 증기압    <dbl> 2.1, 2.2, 2.2, 2.4, 2.2, 2.4, 2.3, 2.3, 2.3, 2.5, 2.6~
$ 이슬점    <dbl> -14.0, -13.4, -13.5, -12.4, -13.3, -12.2, -12.6, -13.~
$ 기압      <dbl> 1020.3, 1020.3, 1020.7, 1020.6, 1020.4, 1020.7, 1021~
$ 해면기압  <dbl> 1024.0, 1024.0, 1024.5, 1024.4, 1024.2, 1024.5, 1025.1~
$ 일조      <dbl> NA, NA, NA, NA, NA, NA, NA, 0.0, 0.0, 0.6, 0.9, 0.5,~
$ 일사      <dbl> 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.0~
$ 전운량    <int> 6, NA, 6, 6, 6, 6, 6, 3, 6, 6, 7, 8, 8, 8, 8, 6, 6, 8~
$ 시정      <dbl> 1500.000, 1750.000, 2000.000, 2000.000, 2000.000, 20~
$ 지면온도  <dbl> -4.4, -4.6, -4.7, -5.0, -5.0, -5.3, -5.7, -5.7, -5.4, ~
$ floating  <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA~
$ warehouse <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 12, 150, 195, 270, 254, 18~
$ dangjin   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 8, 227, 263, 356, 333, 225~
$ ulsan     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
$ hour      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,~

skim(rdata)

Table 1: Data summary

Name	rdata
Number of rows	55392
Number of columns	20
Key	NULL
_______________________
Column type frequency:
numeric	19
POSIXct	1
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
기온	0	1.00	12.21	10.36	-19.3	3.60	12.60	21.10	36.80	▁▅▇▇▂
강수	0	1.00	0.12	1.05	0.0	0.00	0.00	0.00	102.70	▇▁▁▁▁
풍속	0	1.00	1.96	1.56	0.0	0.70	1.60	2.90	11.70	▇▃▁▁▁
풍향	0	1.00	164.29	129.50	0.0	20.00	180.00	290.00	360.00	▇▁▃▂▆
습도	0	1.00	75.34	20.14	10.0	61.00	79.00	94.00	100.00	▁▂▃▅▇
증기압	82	1.00	12.93	8.92	1.0	5.40	10.20	19.20	42.90	▇▅▃▂▁
이슬점	85	1.00	7.37	11.04	-22.4	-1.70	7.30	16.90	30.20	▁▆▇▇▅
기압	79	1.00	1014.05	8.36	983.6	1007.40	1014.30	1020.60	1036.30	▁▃▇▇▂
해면기압	78	1.00	1017.37	8.48	986.6	1010.60	1017.60	1024.00	1039.70	▁▃▇▇▂
일조	25453	0.54	0.52	0.45	0.0	0.00	0.60	1.00	1.00	▇▁▁▁▇
일사	0	1.00	0.54	0.83	0.0	0.00	0.01	0.89	4.85	▇▂▁▁▁
전운량	12427	0.78	5.23	3.84	0.0	1.00	6.00	9.00	10.00	▇▂▃▅▇
시정	0	1.00	1754.20	1024.17	3.0	1100.00	1800.00	2029.00	6454.00	▅▇▁▁▁
지면온도	0	1.00	12.21	10.36	-19.3	3.60	12.60	21.10	36.80	▁▅▇▇▂
floating	28344	0.49	122.42	192.34	0.0	0.00	0.00	192.25	753.00	▇▁▁▁▁
warehouse	2040	0.96	95.81	150.72	0.0	0.00	0.00	152.00	595.00	▇▁▁▁▁
dangjin	2040	0.96	140.11	221.67	0.0	0.00	0.00	227.00	881.00	▇▁▁▁▁
ulsan	3528	0.94	66.28	104.17	0.0	0.00	0.00	104.00	448.00	▇▁▁▁▁
hour	0	1.00	11.50	6.92	0.0	5.75	11.50	17.25	23.00	▇▇▆▇▇

Variable type: POSIXct

skim_variable	n_missing	complete_rate	min	max	median	n_unique
time	0	1	2015-01-01 01:00:00	2021-04-27	2018-02-28 00:30:00	55392

Data preprocessing

rdata %>% 
    select(-hour) %>% 
    mutate(time = ymd_hms(time)) %>% 
    filter(between(time, ymd('2018-03-01'), ymd('2021-01-31'))) -> rdata
rdata %>% glimpse()  

Rows: 25,609
Columns: 19
$ time      <dttm> 2018-03-01 00:00:00, 2018-03-01 01:00:00, 2018-03~
$ 기온      <dbl> 3.1, 2.8, 2.6, 2.0, 2.2, 4.1, 3.5, 2.2, 1.0, 0.3, 0.~
$ 강수      <dbl> 0.50, 0.00, 0.00, 0.00, 0.00, 0.00, 1.80, 0.00, 0.00~
$ 풍속      <dbl> 3.6, 0.7, 3.2, 1.9, 2.1, 4.4, 7.9, 6.4, 7.7, 8.9, 7.~
$ 풍향      <dbl> 340, 140, 320, 230, 180, 270, 320, 290, 320, 320, 32~
$ 습도      <dbl> 96, 97, 95, 97, 97, 97, 93, 86, 82, 71, 63, 58, 60, ~
$ 증기압    <dbl> 7.3, 7.2, 7.0, 6.8, 6.9, 7.9, 7.3, 6.1, 5.4, 4.4, 4.0~
$ 이슬점    <dbl> 2.5, 2.3, 1.8, 1.5, 1.7, 3.6, 2.4, 0.0, -1.7, -4.3, -~
$ 기압      <dbl> 1001.3, 1001.9, 1002.6, 1002.8, 1003.0, 1001.8, 1002~
$ 해면기압  <dbl> 1004.9, 1005.5, 1006.2, 1006.4, 1006.6, 1005.4, 1006.3~
$ 일조      <dbl> NA, NA, NA, NA, NA, NA, NA, NA, 0.0, 0.8, 1.0, 1.0, ~
$ 일사      <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.05~
$ 전운량    <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
$ 시정      <dbl> 922, 4315, 2601, 1717, 1957, 571, 57, 436, 593, 593,~
$ 지면온도  <dbl> 3.1, 2.8, 2.6, 2.0, 2.2, 4.1, 3.5, 2.2, 1.0, 0.3, 0.6,~
$ floating  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 36, 313, 532, 607, 614,~
$ warehouse <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 33, 209, 296, 315, 474,~
$ dangjin   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 37, 318, 490, 550, 727,~
$ ulsan     <int> 0, 0, 0, 0, 0, 0, 0, 0, 4, 35, 71, 82, 334, 372, 3~

rdata %>% 
  summarise(across(.fns = ~sum(is.na(.))/length(.)))

  time 기온 강수 풍속 풍향 습도      증기압     이슬점        기압
1    0    0    0    0    0    0 0.001835292 0.00191339 0.001718146
     해면기압      일조 일사    전운량 시정 지면온도 floating
1 0.001679097 0.4545668    0 0.1552579    0        0        0
  warehouse dangjin ulsan
1         0       0     0

Univariate timeseries analysis

울산 지역의 전력 발전량 데이터만 활용할 것이기 때문에 날짜 변수를 제외한 나머지 변수는 제거했다.

ulsan <- rdata %>% 
  select(-c(dangjin, warehouse, floating)) %>% 
  select(time, ulsan) %>% 
  rename(date = time, value = ulsan)

Time series visualization

tidymodels의 경우 train/test 분리를 위해서 initial_split()을 활용했는데 시계열 데이터의 경우 특정 날짜를 기준으로 잘라야하기 때문에 Modeltime 패키지에 내장되어있는 initial_time_split()을 이용한다. 특정 날짜를 기준으로 자르고 싶을 경우 timeseries_split()을 이용할 수도 있다.

• tk_time_series_cv_plan() : split된 object를 데이터프레임으로 전환

• plot_time_series_cv_plan() : sampling된 데이터를 이용해서 시계열 그래프 생성

initial_time_split(data = ulsan, prop = 0.9) %>% 
  tk_time_series_cv_plan() %>% 
  plot_time_series_cv_plan(date, value,
                           .interact = FALSE, 
                           .title = "Partition Train / Test")

Split train/test

months <- 1

total_months <- lubridate::interval(base::min(ulsan$date),
                                    base::max(ulsan$date)) %/%  
                                    base::months(1)


prop <- (total_months - months) / total_months

splits <- rsample::initial_time_split(ulsan, prop = prop)


splits %>%
  timetk::tk_time_series_cv_plan() %>%  
  timetk::plot_time_series_cv_plan(date, value) 

Model fitting

모델 fitting은 tidymodels 패키지의 방식과 동일하다. 현재는 default 세팅으로 모델을 fitting했지만 튜닝 파라미터가 있을 경우 이전에 tidymodels에서 했던 방식 그대로 grid_latin_hypercube(), grid_random, bayes tuning 등을 이용해서 최적의 파라미터를 찾고 모델을 fitting 해야한다.

# Exponential smoothing 
model_fit_ets <- modeltime::exp_smoothing() %>%
    parsnip::set_engine(engine = "ets") %>%
    parsnip::fit(value ~ date, data = training(splits))


# ARIMA 
model_fit_arima <- modeltime::arima_reg() %>%
    parsnip::set_engine("auto_arima") %>%
    parsnip::fit(
        value ~ date, 
        data = training(splits))


# ARIMA Boost
model_fit_arima_boost <- modeltime::arima_boost() %>%
    parsnip::set_engine("auto_arima_xgboost") %>%
    parsnip::fit(
        value ~ date + as.numeric(date) + month(date), 
        data = training(splits))

# Prophet
model_fit_prophet <- modeltime::prophet_reg() %>%
    parsnip::set_engine("prophet") %>%
    parsnip::fit(
        value ~ date, 
        data = training(splits))


# Prophet Boost
model_fit_prophet_boost <- modeltime::prophet_boost() %>%
    parsnip::set_engine("prophet_xgboost") %>%
    parsnip::fit(
        value ~ date + as.numeric(date) + month(date), 
        data = training(splits))

Model time table

modeltime_table에 fitting한 모델을 추가한다. modeltime_table은 이전에 서술했다시피 각 모델이 재대로 적합되었는지 확인하고, 이후 예측 워크플로우를 위해서 modeltime_table 구조를 이용하므로 model table을 에러 없이 세팅하는 것이 중요하다.

model_tbl <- modeltime::modeltime_table(
    model_fit_ets,
    model_fit_arima,
    model_fit_arima_boost,
    model_fit_prophet,
    model_fit_prophet_boost)

model_tbl

# Modeltime Table
# A tibble: 5 x 3
  .model_id .model   .model_desc                              
      <int> <list>   <chr>                                    
1         1 <fit[+]> ETS(A,N,A)                               
2         2 <fit[+]> ARIMA(5,0,0)(2,1,0)[24]                  
3         3 <fit[+]> ARIMA(2,0,1)(2,1,0)[24] W/ XGBOOST ERRORS
4         4 <fit[+]> PROPHET                                  
5         5 <fit[+]> PROPHET W/ XGBOOST ERRORS                

calibration

이전에 만든 modeltime_table을 test 데이터에 적합시켜서 보정을 하는 단계이다.

calibration_tbl <- model_tbl %>%
    modeltime::modeltime_calibrate(testing(splits))  

calibration_tbl %>%
    modeltime::modeltime_accuracy() %>%   
    flextable::flextable() %>% 
    flextable::bold(part = "header") %>% 
    flextable::bg(bg = "#D3D3D3", part = "header") %>% 
    flextable::autofit()

visualization the forecast test

calibration_tbl %>%
    modeltime::modeltime_forecast(new_data = testing(splits), 
                                  actual_data = ulsan,
                                  conf_interval = 0.90) %>%
    modeltime::plot_modeltime_forecast(.legend_show = TRUE, 
                                       .legend_max_width = 20)

Refit

train/test를 합한 original 데이터를 이용해서 refitting을 진행하는 단계이다.

refit_tbl <- calibration_tbl %>%
    modeltime::modeltime_refit(data = ulsan)

forecast

refitting된 모델을 이용해서 지정한 time interval에 대해 forecast를 수행하는 단계이다.

forecast_tbl <- refit_tbl %>%
    modeltime::modeltime_forecast(
        h = "1 month",
        actual_data = ulsan,
        conf_interval = 0.90
    ) 

forecast_tbl %>%
    modeltime::plot_modeltime_forecast(.interactive = TRUE,
                                       .legend_max_width = 20)

Aggregate model

Accuracy 향상을 위해서 적합시킨 5 가지 모델을 평균내서 최종 모델을 산출한다. Modeltime.ensemble을 이용하면 더 세련된 앙상블 기법을 이용할 수 있다.

mean_forecast_tbl <- forecast_tbl %>%
    dplyr::filter(.key != "actual") %>%
    dplyr::group_by(.key, .index) %>%
    dplyr::summarise(across(.value:.conf_hi, mean)) %>%
    dplyr::mutate(
        .model_id   = 6,
        .model_desc = "AVERAGE OF MODELS"
    )


# Visualize aggregate model 
forecast_tbl %>%
    dplyr::filter(.key == "actual") %>%
    dplyr::bind_rows(mean_forecast_tbl) %>%
    modeltime::plot_modeltime_forecast()  

참고 자료

H2O 관련

https://www.r-bloggers.com/2021/03/introducing-modeltime-h2o-automatic-forecasting-with-h2o-automl/

gluonTS 관련

https://cran.r-project.org/web/packages/modeltime.gluonts/vignettes/getting-started.html#references

Modeltime 관련

https://cran.r-project.org/web/packages/modeltime/vignettes/getting-started-with-modeltime.html

https://www.adam-d-mckinnon.com/posts/2020-10-10-forecastpeopleanalytics/